home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-02
/
stat.zip
/
STAT.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1990-08-08
|
45KB
|
1,590 lines
{--------------------------------------------------------------------------}
{ Norton Statistical Library }
{ }
{ Version 1.00 }
{ }
{ }
{ Copyright 1990 Norton Associcates }
{ All Rights Reserved }
{ Restricted by License }
{--------------------------------------------------------------------------}
{--------------------------------}
{ Unit: Stat }
{--------------------------------}
{$S-,R-,V-,D-,A+,B+,N+,E-,I-}
UNIT
stat;
INTERFACE
CONST
error_value = -99.9; { if any errors }
maxsize = 65520; { max segment size }
maxsingle = maxsize DIV (SIZEOF(SINGLE)); { max element size single array }
maxdouble = maxsize DIV (SIZEOF(DOUBLE)); { max element size double array }
maxlongint = maxsize DIV (SIZEOF(LONGINT)); { max element size longint array }
v1 : INTEGER = 1; { constants for uniform1 }
v2 : INTEGER = 1000; { constants for uniform1 }
v3 : INTEGER = 30000; { constants for uniform1 }
maxorder = 10;
maxmatrix = maxorder * 2 - 1;
uno = 1;
dos = 2;
zero = 0.0;
one = 1.0;
two = 2.0;
c2 = 2.0;
c3 = -3.0;
c4 = 4.0;
c5 = 5.0;
c6 = 6.0;
c11 = 11.0;
c12 = 12.0;
c17 = 17.0;
i_4 = 1.0/c4;
i_16 = 1.0/16.0;
i_35 = 1.0/35.0;
TYPE
single_array_type = SINGLE;
single_array_dummy = ARRAY[1..maxsingle] OF single_array_type;
single_array_pointer = ^single_array_dummy;
double_array_dummy = ARRAY[1..maxdouble] OF DOUBLE;
double_array_pointer = ^double_array_dummy;
longint_array_dummy = ARRAY[1..maxlongint] OF LONGINT;
longint_array_pointer = ^longint_array_dummy;
quartype = ARRAY[1..5] OF SINGLE;
arry_type = ARRAY[1..maxorder,1..maxorder] OF EXTENDED;
PROCEDURE create_single_array( num: WORD; VAR xx : single_array_pointer);
PROCEDURE delete_single_array( num: WORD; VAR xx : single_array_pointer);
PROCEDURE create_longint_array( num: WORD; VAR xx : longint_array_pointer);
PROCEDURE delete_longint_array( num: WORD; VAR xx : longint_array_pointer);
FUNCTION uniform1 : SINGLE;
FUNCTION rndnorm1( mean,standev:EXTENDED) :SINGLE;
FUNCTION rndnorm2( mean,standev:EXTENDED) :SINGLE;
PROCEDURE insert( n : WORD ; VAR a : single_array_pointer) ;
PROCEDURE qsort( n : WORD; VAR a : single_array_pointer) ;
PROCEDURE remove_avg( n : WORD; VAR a : single_array_pointer ; avg : SINGLE);
PROCEDURE means( n :WORD; a : single_array_pointer; VAR xmean,gmean,hmean,rmsmean :SINGLE);
PROCEDURE wxmean( num : WORD; a : single_array_pointer; freq : longint_array_pointer; VAR mean,sd,small,large : SINGLE);
PROCEDURE elem_stat( num:WORD; VAR a:single_array_pointer;
VAR small,large:SINGLE; VAR mean,sd:SINGLE);
PROCEDURE moments( n : WORD;
a : single_array_pointer;
VAR ave : SINGLE;
VAR std : SINGLE;
VAR skew : SINGLE;
VAR beta2 : SINGLE);
PROCEDURE quart( n : WORD; a : single_array_pointer; VAR quart : quartype);
FUNCTION percentile( n : WORD; a : single_array_pointer ;percent : SINGLE) : SINGLE;
FUNCTION standard_error(num : WORD ; sd : SINGLE; ntype:WORD) : SINGLE;
FUNCTION cdf_prob_to_sd(prob :SINGLE) : SINGLE;
FUNCTION cdf_sd_to_prob(sd:SINGLE) : SINGLE;
FUNCTION int_prob_to_sd(prob :SINGLE) : SINGLE;
FUNCTION int_sd_to_prob(sd:SINGLE) : SINGLE;
PROCEDURE corcoef(n: WORD; x,y:single_array_pointer; VAR r:SINGLE);
PROCEDURE autocor(n:WORD; x:single_array_pointer; lag :WORD; VAR auto:single_array_pointer);
PROCEDURE linfit(npts: WORD; x,y,sigmay: single_array_pointer; mode:WORD;VAR a, b, r:SINGLE);
FUNCTION determ(arry :arry_type; norder : INTEGER) : EXTENDED;
PROCEDURE polfit(npts: WORD;
x,y,sdy : single_array_pointer;
nterms,mode : INTEGER;
VAR a : single_array_pointer;
VAR r : SINGLE;
VAR se : SINGLE);
PROCEDURE mulreg(n : WORD; y,x,z : single_array_pointer;
VAR a : single_array_pointer;
VAR r : SINGLE;
VAR se : SINGLE);
PROCEDURE smooth121(n:WORD; VAR y: single_array_pointer);
PROCEDURE smooth14641(n:WORD; VAR y: single_array_pointer);
PROCEDURE smoothcurve(n:WORD; VAR y: single_array_pointer);
FUNCTION movavg(n: WORD; a : single_array_pointer; ma:WORD ; k : WORD) : SINGLE;
{*****************************************************************************}
{*****************************************************************************}
IMPLEMENTATION
{*****************************************************************************}
{*****************************************************************************}
PROCEDURE create_single_array( num: WORD; VAR xx : single_array_pointer);
{ Author : Norton Associates
Purpose: Properly create a single dimensioned heap array
Version: 1.0
Date : 5 May 1990 }
VAR
maxnumber : LONGINT; { proposed size of array in bytes }
BEGIN
maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
IF (num < uno) or (num > maxsingle) THEN
BEGIN
WRITELN('Sorry the single array size is to large to create ');
WRITELN('You wanted = ',num:10,', The max single size = ',maxsingle:10);
HALT;
END;
GETMEM(xx,maxnumber);
END;
PROCEDURE delete_single_array( num: WORD; VAR xx : single_array_pointer);
{ Author : Norton Associates
Purpose: Properly delete a single dimensioned heap array
Version: 1.0
Date : 5 May 1990 }
VAR
maxnumber : LONGINT; { proposed size of array in bytes }
BEGIN
maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
IF maxnumber > maxsize THEN
BEGIN
WRITELN('sorry the single array size is to large to delete ',maxnumber);
HALT;
END;
FREEMEM(xx,maxnumber);
END;
PROCEDURE create_longint_array( num: WORD; VAR xx : longint_array_pointer);
{ Author : Norton Associates
Purpose: Properly create a longint dimensioned heap array
Version: 1.0
Date : 5 May 1990 }
VAR
maxnumber : LONGINT; { proposed size of array in bytes }
BEGIN
maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
IF (num < uno) or (num > maxlongint) THEN
BEGIN
WRITELN('sorry longint array size is to large to create',maxlongint);
HALT;
END;
GETMEM(xx,maxnumber);
END;
PROCEDURE delete_longint_array( num: WORD; VAR xx : longint_array_pointer);
{ Author : Norton Associates
Purpose: Properly delete a longint dimensioned heap array
Version: 1.0
Date : 5 May 1990 }
VAR
maxnumber : LONGINT; { proposed size of array in bytes }
BEGIN
maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
IF maxnumber > maxsize THEN
BEGIN
WRITELN('sorry longint array size is to large to delete',maxnumber);
HALT;
END;
GETMEM(xx,maxnumber);
END;
PROCEDURE remove_avg(n : WORD; VAR a : single_array_pointer ; avg : SINGLE);
{ Author : Norton Associates
Purpose: Remove a constant value from an array
Version: 1.0
Date : 5 May 1990 }
VAR
j : WORD;
BEGIN
IF n > 0 THEN
BEGIN
FOR j := uno TO n DO
a^[j] := a^[j] - avg;
END
ELSE
WRITELN('remavg> n < 1: No data');
END;
PROCEDURE wxmean(num : WORD; a : single_array_pointer; freq : longint_array_pointer; VAR mean,sd,small,large : SINGLE);
{ Author : Norton Associates
Purpose: Calculate a weighted arithemetic mean
Version: 1.0
Date : 5 May 1990 }
VAR
number, { num = number of variables to process }
j : WORD; { loop index }
dum,
sum,sum2 : EXTENDED; { sum of squares }
BEGIN
{ check to see if there is enough data }
IF num < dos THEN
BEGIN
WRITELN('Wxmean> No sample data: n < 2');
small := error_value;
large := error_value;
mean := error_value;
sd := error_value;
EXIT;
END;
{ zero out and initalize some data }
sum := zero;
sum2 := zero;
small := a^[1];
large := a^[1];
{ calculate sum of squares }
FOR j := uno TO num DO
BEGIN
dum := a^[j] * freq^[j];
sum := sum + dum;
sum2 := sum2 + SQR(dum);
number := number + freq^[j];
IF a^[j] > large THEN
large := a^[j]
ELSE IF a^[j] < small THEN
small := a^[j];
END;
{ calculate mean of sample population }
mean := sum / number;
{ check to see if standard deviation can be calculated }
IF number - uno < uno THEN
BEGIN
sd := error_value;
WRITELN('wxmean> no standard deviation calculated');
EXIT;
END
ELSE IF (sum2 - (number * mean * mean) ) <= zero THEN
sd := zero
ELSE
sd := SQRT((sum2 - (number * mean * mean) ) / (number-1));
END;
PROCEDURE means(n :WORD; a : single_array_pointer; VAR xmean,gmean,hmean,rmsmean:SINGLE);
{ Author : Norton Associates
Purpose: Calculate the arithemetic, geometric, root mean square and
harmonic mean
Version: 1.0
Date : 5 May 1990 }
VAR
xsum,gsum,hsum,rmssum : EXTENDED;
gboolean,hboolean : BOOLEAN;
j : WORD;
BEGIN
{ clear and set up the data }
xsum := zero; { arithmetic mean }
hsum := zero; { harmonic mean }
gsum := zero; { geometric mean }
rmssum := zero; { roor mean square mena }
gboolean := TRUE;
hboolean := TRUE;
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('means> no sample data: n < 2');
EXIT;
END;
{ sort out the data }
FOR j := uno TO n DO
BEGIN
xsum := xsum + a^[j];
IF a^[j] > zero THEN
gsum := gsum + LN(a^[j])
ELSE
gboolean := FALSE; { bad data flag }
IF a^[j] > zero THEN
hsum := hsum + one / a^[j]
ELSE
hboolean := FALSE; { bad data flag }
rmssum := rmssum + SQR(a^[j]);
END;
IF gboolean = FALSE THEN
BEGIN
WRITELN('means> Gmean can not be calculated: Hmean = -99.9');
gmean := error_value
END
ELSE
gmean := EXP(gsum/n);
IF hboolean = FALSE THEN
BEGIN
WRITELN('means> Hmean can not be calculated: Hmean = -99.9');
hmean := error_value
END
ELSE
gmean := EXP(gsum/n);
xmean := xsum / n;
rmsmean := SQRT(rmssum / n);
END;
FUNCTION cdf_prob_to_sd(prob :SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the standard deviation from probability
Version: 1.0
Date : 5 May 1990
PROB=PROBABILITY (0-1) INPUT
SD=NO. OF STANDARD DEVIATIONS,OUTPUT }
CONST
c0 : extended = 2.515517;
c1 : extended = 0.802853;
c2 : extended = 0.010328;
d1 : extended = 1.432788;
d2 : extended = 0.189269;
d3 : extended = 0.001308;
VAR
dn,up,xp : EXTENDED;
t1,t2,t3 : EXTENDED;
BEGIN
{ check to see if the probability ( prob ) is in range }
IF (prob < zero) OR (prob > one) THEN
BEGIN
WRITELN('prob_to_sd> input probability out of range (0-1)');
EXIT;
END;
IF prob > 0.5 THEN
xp := one - prob
ELSE
xp := prob;
t1 := SQRT(LN(one/SQR(xp)));
t2 := t1 * t1;
t3 := t2 * t1;
up := c0 + c1*t1 + c2*t2;
dn := one + d1*t1 + d2*t2 + d3*t3;
xp := t1 - (up/dn);
IF prob <= 0.5 THEN
cdf_prob_to_sd := xp * (-one);
cdf_prob_to_sd := xp;
END;
FUNCTION cdf_sd_to_prob(sd:SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the standard deviation from probability
Version: 1.0
Date : 5 May 1990
sd = (mean-a)/standard deviation
prob = percentile of a between 0 and 100%)
}
CONST
p : extended = 0.231641900;
b1 : extended = 0.319381530;
b2 : extended = -0.356563782;
b3 : extended = 1.781477937;
b4 : extended = -1.821255978;
b5 : extended = 1.330274429;
VAR
x,y1,y2,y3,y4,y5 : EXTENDED;
t,z,r : EXTENDED;
BEGIN
{ check to see if the standard deviation ( sd ) is in range }
IF sd < zero THEN
BEGIN
WRITELN('cdf_sd_prob> input standard deviation out of range: sd >=0.0');
EXIT;
END;
x := sd;
r := EXP(-(x*x)/two)/2.5066282746;
{ r = frequency }
z := x;
y1 := one/(one+(p*ABS(x)));
y2 := y1 * y1;
y3 := y2 * y1;
y4 := y3 * y1;
y5 := y4 * y1;
t := one - r*(b1*y1+ b2*y2 + b3*y3 + b4*y4 + b5*y5);
IF z > zero THEN
cdf_sd_to_prob := t
ELSE
cdf_sd_to_prob := one - t;
{ prob = [0% to 100%]}
END;
FUNCTION int_prob_to_sd(prob :SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the probability from the standard deviation
Version: 1.0
Date : 5 May 1990 }
VAR
dummy : SINGLE;
BEGIN
{ check to see if the probability is in range }
IF (prob < zero) OR (prob > one) THEN
BEGIN
WRITELN('int_prob_to_sd> input probability out of range: prob (0-1)');
EXIT;
END;
dummy := 0.50 + (0.50 * prob) ;
int_prob_to_sd := cdf_prob_to_sd(dummy);
END;
FUNCTION int_sd_to_prob(sd:SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the probability from the standard deviation
Version: 1.0
Date : 5 May 1990 }
VAR
dummy: SINGLE;
BEGIN
{ check to see if the standard deviation ( sd ) is in range }
IF sd <= zero THEN
BEGIN
WRITELN('int_sd_to_prob> negative or zero standard deviation not allowed');
EXIT;
END;
dummy := cdf_sd_to_prob(sd);
int_sd_to_prob := dummy - (0.50 - (dummy - 0.50) )
END;
PROCEDURE linfit(npts : WORD; x,y,sigmay: single_array_pointer;
mode : WORD; VAR a, b, r : SINGLE);
{ Author : Norton Associates
Purpose: Calculate a linear least squares fit for x and y pairs of data
Version: 1.0
Date : 5 May 1990 }
VAR
sum,sumx,sumy,sumx2,sumy2,sumxy : EXTENDED;
x1,y1,weight,delta : EXTENDED;
i : WORD;
{ x,y pairs of data
sigmay standard deviation of y
npts = number of pairs of as
mode = 0 = no weighting
mode = 1 = instrumental weighting (sigmay array must exist)
a = intercept
b = slope
r = correlation coefficient }
BEGIN
{ zero out items for use }
sum := zero;
sumx := zero;
sumy := zero;
sumxy := zero;
sumx2 := zero;
sumy2 := zero;
{ sum over npts pairs of points }
{ check to see if there is enough data }
IF npts < dos THEN
BEGIN
WRITELN('linfit> number of data points < 2');
EXIT;
END;
{ calculate sum of square for the data }
FOR i := uno TO npts DO
BEGIN
x1 := x^[i];
y1 := y^[i];
IF (mode = 0) OR (y1 = zero) THEN
weight := one
ELSE
weight := one / SQR(sigmay^[i]);
sum := sum + weight;
sumx := sumx + weight * x1;
sumy := sumy + weight * y1;
sumx2 := sumx2 + weight * x1 * x1;
sumy2 := sumy2 + weight * y1 * y1;
sumxy := sumxy + weight * x1 * y1;
END;
{ calculate final answers }
delta := sum * sumx2 - sumx * sumx;
a := (sumx2*sumy - sumx*sumxy) / delta;
b := (sumxy*sum - sumx*sumy) / delta;
r := (sum*sumxy - sumx*sumy) /
SQRT(delta * (sum * sumy2 - sumy * sumy));
END;
FUNCTION standard_error(num : WORD ; sd : SINGLE; ntype : WORD) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the standard error based upon a normal distribution
Version: 1.0
Date : 5 May 1990 }
BEGIN
{ check to see if there is enough data }
IF num < dos THEN
BEGIN
WRITELN('Standard_error> number of samples < 2');
EXIT;
END;
{ determine the type and calculate the standard error }
CASE ntype OF
1 : standard_error := sd / SQRT(num);
2 : standard_error := sd * SQRT(num);
3 : standard_error := 1.2533 * sd / SQRT(num);
4 : standard_error := sd / SQRT(two * num);
5 : standard_error := 0.6028 * sd / SQRT(num);
6 : standard_error := sd * SQRT(one + (two * SQR(sd))) / SQRT(two * num);
7 : standard_error := (one - SQR(sd)) / SQRT(num);
END;
END;
PROCEDURE qsort( n : WORD ; { number of as }
VAR a : single_array_pointer) ; { number of as in a }
{ Author : Norton Associates
Purpose: Sort data based upon quick sort algorithm
Version: 1.0
Date : 5 May 1990 }
CONST
brkeven = 13; { more than 10 as use qsort otherwise use insert sort}
maxstack = 20; { enough stack space for 1_000_000 random numbers }
VAR
low,high,ti,median,jd,
bi,id,index,pt,j,i : WORD;
stack : 0..maxstack;
top,bottom : ARRAY[1..maxstack] OF WORD;
pivot : SINGLE; { change to match singlearray }
PROCEDURE swap_position( median,low : WORD);
VAR
t : SINGLE;
BEGIN
t := a^[median];
a^[median] := a^[low];
a^[low] := t;
END;
BEGIN { main program }
{ non-recursive qsort, must use stack : 10% performance increase }
{ initialize the qsort stack }
stack := uno;
top[stack] := uno;
bottom[stack] := n;
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('qsort> no sample data');
EXIT;
END;
REPEAT
{ pop information off the stack }
low := top[stack];
high := bottom[stack];
{ use of insert sort for small sublists : 16% performance increase }
{ use qsort for more than brkeven or use insert sort for small sublists }
WHILE (high - low) > brkeven DO
BEGIN
ti := low;
bi := high;
{ median of three rule : 10-100% performance increase, almost no worst case }
median := (low + high) DIV 2; { initial guess }
IF a^[low] > a^[median] THEN
swap_position(median,low);
IF a^[high] < a^[median] THEN
BEGIN
swap_position(high,median);
IF a^[low] > a^[median] THEN
swap_position(low,median);
END;
pivot := a^[median]; { pivot = median of three }
{ the old qsort partition algorithm, median of three optimized}
REPEAT
REPEAT dec(bi); UNTIL a^[bi] <= pivot;
REPEAT inc(ti); UNTIL a^[ti] >= pivot;
IF ti <= bi THEN
swap_position(bi,ti);
UNTIL ti > bi ;
{ larger of the two stacks saved first: 10% performance increase}
IF (bi - low) > (high - ti) THEN
BEGIN
top[stack] := low;
bottom[stack] := bi;
low := ti;
END
ELSE
BEGIN
top[stack] := ti;
bottom[stack] := high;
high := bi;
END;
{ wrap up non-recursive part of qsort }
inc(stack);
IF stack > maxstack THEN
BEGIN
WRITELN('qsort> stack space exceeded');
HALT;
END;
END;
{ use insertion sort for small sublists : 16% performance increase}
IF (high - low) > 0 THEN
BEGIN
id := low + uno;
FOR i := id TO high DO
BEGIN
pivot := a^[i];
IF a^[low] > pivot THEN
BEGIN
FOR j := i DOWNTO id DO
a^[j] := a^[j-1];
a^[low] := pivot;
END
ELSE
BEGIN
j := i;
WHILE a^[j-1] > pivot DO
BEGIN
a^[j] := a^[j-1];
dec(j);
END;
a^[j] := pivot;
END
END;
END;
dec(stack);
UNTIL stack = 0
END;
PROCEDURE insert(n : WORD ; { number of as }
VAR a : single_array_pointer) ; { number of as in a }
{ Author : Norton Associates
Purpose: Sort data based upon insertion sort algorithm
Version: 1.0
Date : 5 May 1990 }
VAR
i,j : WORD;
pivot : SINGLE;
BEGIN
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('insert> no sample data');
EXIT;
END;
FOR i := 2 TO n DO
BEGIN
pivot := a^[i];
IF a^[1] > pivot THEN
BEGIN
FOR j := i DOWNTO 2 DO
a^[j] := a^[j-1];
a^[1] := pivot;
END
ELSE
BEGIN
j := i;
WHILE a^[j-1] > pivot DO
BEGIN
a^[j] := a^[j-1];
dec(j);
END;
a^[j] := pivot;
END
END;
END;
FUNCTION uniform1 : SINGLE;
{ Author : Norton Associates
Purpose: Sort data based upon insertion sort algorithm
Version: 1.0
Date : March 1987 Byte Magazine }
VAR
temp : SINGLE ;
BEGIN
v1 := 171 * ( v1 MOD 177) - 2 * ( v1 DIV 177);
IF v1 < 0 THEN
v1 := v1 + 30269;
v2 := 172 * (v2 MOD 176) - 35 * ( v2 DIV 176);
IF v2 < 0 THEN
v2 := v2 + 30307;
v3 := 170 * (v3 MOD 178) - 63 * ( v3 MOD 178);
IF v3 < 0 THEN
v3 := v3 + 30323;
temp := (v1 / 30269.0) + (v2 / 30307.0) + (v3 / 30323.0);
uniform1 := temp - TRUNC(temp);
END;
FUNCTION rndnorm1(mean , standev : EXTENDED) :SINGLE;
{ Author : Norton Associates
Purpose: Calculate random normal number
Version: 1.0
Date : 5 May 1990 }
VAR
randoma,randomb,randomc,radius2,deviate : EXTENDED;
BEGIN
REPEAT
randoma := two * RANDOM - one;
randomb := two * RANDOM - one;
radius2 := SQR(randoma) + SQR(randomb);
UNTIL radius2 < one;
IF RANDOM(2) = uno THEN
randomc := randomb
ELSE
randomc := randoma;
deviate := randomc * SQRT(( - two * LN(radius2)) / radius2);
rndnorm1 := mean + deviate * standev;
END;
FUNCTION rndnorm2( mean , standev : EXTENDED) :SINGLE;
{ Author : Norton Associates
Purpose: Calculate random normal number ( approximation to rndnom1 )
Version: 1.0
Date : 5 May 1990 }
VAR
sum : SINGLE;
j : WORD;
BEGIN
sum := zero;
FOR j := uno TO 12 DO
sum := sum + RANDOM;
rndnorm2 := mean + (sum - 6.0) * standev;
END;
FUNCTION power( number, exponent : SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Raise a number to a number
Version: 1.0
Date : 5 May 1990 }
BEGIN
IF (number <> zero) AND (exponent <> zero) THEN
power := EXP( exponent * LN(number) )
ELSE
WRITELN('power> can not take power of ',number:10,exponent:10);
END;
PROCEDURE elem_stat(num:WORD; VAR a:single_array_pointer;
VAR small,large:SINGLE; VAR mean,sd:SINGLE);
{ Author : Norton Associates
Purpose: Calculate mean, standard deviation the largest and smallest value
Version: 1.0
Date : 5 May 1990 }
VAR
j : WORD; { loop index }
sum,sum2 : EXTENDED; { sum of squares }
BEGIN
{ check to see if there is data }
IF num < uno THEN
BEGIN
WRITELN('elem_stat> no data points');
small := error_value;
large := error_value;
mean := error_value;
sd := error_value;
EXIT;
END;
{ zero out data }
sum := zero;
sum2 := zero;
small := a^[1];
large := a^[1];
{ calculate sum of squares }
FOR j := uno TO num DO
BEGIN
sum := sum + a^[j];
sum2 := sum2 + SQR(a^[j]);
IF a^[j] > large THEN
large := a^[j]
ELSE IF a^[j] < small THEN
small := a^[j];
END;
{ calculate mean and standard deviation }
mean := sum / num;
IF num - uno < uno THEN
BEGIN
sd := error_value;
WRITELN('elem_stat> number of samples = 1 no standard deviation');
EXIT;
END
ELSE IF (sum2 - (num * mean * mean) ) <= zero THEN
sd := zero
ELSE
sd := SQRT((sum2 - (num * mean * mean) ) / (num-1));
END;
PROCEDURE moments( n : WORD;
a : single_array_pointer;
VAR ave : SINGLE;
VAR std : SINGLE;
VAR skew : SINGLE;
VAR beta2 : SINGLE);
{ Author : Norton Associates
Purpose: Calculate the first 4 moments of the sample population
Version: 1.0
Date : 5 May 1990 }
{ COMPUTE STATISTICAL MOMENTS }
VAR
d,sum,sum2,sum3,sum4,z2,z3,z4,s2,s3,s4 : EXTENDED;
t : SINGLE;
i : WORD;
BEGIN
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('moments> no sample data');
EXIT;
END;
{...COMPUTE MEAN }
t := n;
sum := zero;
FOR i := uno TO n DO
sum := sum + a^[i];
ave := sum / t;
{...COMPUTE OTHER MOMENTS }
sum2 := zero;
sum3 := zero;
sum4 := zero;
FOR i := uno TO n DO
BEGIN
d := a^[i] - ave;
z2 := d * d;
z3 := d * z2;
z4 := d * z3;
sum2 := sum2 + z2;
sum3 := sum3 + z3;
sum4 := sum4 + z4;
END;
{...COMPUTE VARIANCE,STANDARD DEVIATION }
s2 := sum2 /(t-one);
std := SQRT(s2);
{...COMPUTE COEFFICIENT OF SKEWNESS }
IF s2 <> zero THEN
BEGIN
s3 := sum3 / t;
skew := s3 / power(s2,1.5);
{...COMPUTE COEFFICIENT OF KURTOSIS }
s4 := sum4 / t;
beta2 := s4 / (s2*s2);
END
ELSE
BEGIN
WRITELN('moments> no data for beta2 and skewness');
beta2 := error_value;
skew := error_value;
END;
END;
PROCEDURE quart( n : WORD;
a : single_array_pointer;
VAR quart : quartype);
{ Author : Norton Associates
Purpose: Calculate the three quartiles and the smallest and largest
of sample population
Version: 1.0
Date : 5 May 1990 }
CONST
p : ARRAY[1..3] OF SINGLE = (0.25,0.50,0.75);
VAR
r,add,dn,dr :SINGLE;
nt,nb,j : WORD;
BEGIN
{ check to see if there is data }
IF n < dos THEN
BEGIN
WRITELN('quart> no sample data');
EXIT;
END;
qsort(n,a);
quart[1] := a^[1]; { smallest }
quart[5] := a^[n]; { largest }
FOR j := uno TO 3 DO
BEGIN
r := p[j] * n;
nt := ROUND(r + 0.999);
nb := TRUNC(r);
add := zero;
IF (nb = 0) OR (nt = 0) THEN
quart[j+1] := a^[1]
ELSE
BEGIN
IF nb <> nt THEN
BEGIN
dn := r - nb;
dr := a^[nt] - a^[nb];
add := dn * dr;
END;
quart[j+1] := a^[nb] + add;
END;
END;
END;
FUNCTION percentile(n : WORD; a : single_array_pointer ; percent : SINGLE) : SINGLE;
{ Author : Norton Associates
Purpose: Calculate the value in the sample population based upon
percentile provided
Version: 1.0
Date : 5 May 1990 }
VAR
r,add,dn,dr :SINGLE;
nt,nb,j : WORD;
BEGIN
{ check to see if percentile is in range }
IF percent < zero THEN
BEGIN
WRITELN('percentile> error percentile must be >= 0.0');
EXIT;
END;
{ any data }
IF n < one THEN
BEGIN
WRITELN('percentile> no sample data');
EXIT;
END;
percent := percent / 100.0;
qsort(n,a);
r := percent * n;
nt := ROUND(r + 0.999);
nb := TRUNC(r);
add := zero;
IF (nb = 0) OR (nt = 0) THEN
percentile := a^[1]
ELSE
BEGIN
IF nb <> nt THEN
BEGIN
dn := r - nb;
dr := a^[nt] - a^[nb];
add := dn * dr;
END;
percentile := a^[nb] + add;
END;
END;
FUNCTION determ(arry :arry_type; norder : INTEGER) : EXTENDED;
{ Author : Norton Associates
Purpose: Solve simultaneous equations
Version: 1.0
Date : 5 May 1990 }
VAR
det,save : EXTENDED;
j,k,k1,
l,i : WORD;
BEGIN
det := one;
FOR k := uno TO norder DO
BEGIN
IF arry[k,k] = zero THEN
BEGIN
FOR j := k TO norder DO
BEGIN
IF arry[k,j] = zero THEN
BEGIN
WRITELN('determinate> determinate contains a zero');
determ := zero;
EXIT;
END;
END;
FOR i := k TO norder DO
BEGIN
save := arry[i,j];
arry[i,j] := arry[i,k];
arry[i,k] := save;
END;
det := - det;
END;
det := det * arry[k,k];
IF k - norder < 0 THEN
BEGIN
k1 := k + uno;
FOR i := k1 TO norder DO
FOR j := k1 TO norder DO
arry[i,j] := arry[i,j] - arry[i,k] * arry[k,j] / arry[k,k];
END;
END;
determ := det;
END;
PROCEDURE polfit(npts : WORD;
x,y,sdy : single_array_pointer;
nterms,mode : INTEGER;
VAR a : single_array_pointer;
VAR r : SINGLE;
VAR se : SINGLE);
{ Author : Norton Associates
Purpose: Determine from least squares the regression
and correlation coefficients
Version: 1.0
Date : 5 May 1990 }
VAR
wt,xterm,yterm,
delta,srs,sum_y,sum_y2,
ycalc,x_power,residual : EXTENDED;
sumx : ARRAY[1..maxmatrix] OF EXTENDED;
sumy : ARRAY[1..maxorder] OF EXTENDED;
arry : arry_type;
n,j,nmax,k,l : WORD;
xj,yj : SINGLE;
sigmay2,se2,dummy : EXTENDED;
PROCEDURE bad_data(nn : WORD);
VAR
j : WORD;
BEGIN
r := zero;
se := zero;
FOR j := uno TO nterms DO
a^[j] := zero;
WRITELN('polfit> Determinate or correlation approximately 0.0 ',nn);
END;
BEGIN
{ check to see if there is data }
IF npts < dos THEN
BEGIN
WRITELN('polfit> number of data points < 2');
EXIT;
END;
{ check the number of terms versus maximum order }
IF nterms > maxorder THEN
BEGIN
WRITELN('polfit> maxorder for polfit is ',maxorder:10);
EXIT;
END;
{ check to see if number of points is more than number of terms }
IF nterms >= npts THEN
BEGIN
WRITELN('polfit> number of terms must be less than number of points');
EXIT;
END;
{ zero out data }
r := zero;
FOR j := uno TO nterms DO
a^[j] := zero;
nmax := 2 * nterms - uno;
FOR n := uno TO nmax DO
sumx[n] := zero;
{ calculate sum of squares for problem }
FOR j := uno TO nterms DO
sumy[j] := zero;
FOR j := uno TO npts DO
BEGIN
xj := x^[j];
yj := y^[j];
IF (mode = 0) OR (yj = zero) THEN
wt := one
ELSE
wt := one/(SQR(sdy^[j]));
xterm := wt;
FOR n := uno TO nmax DO
BEGIN
sumx[n] := sumx[n] + xterm;
xterm := xterm * xj;
END;
yterm := wt * yj;
FOR n := uno TO nterms DO
BEGIN
sumy[n] := sumy[n] + yterm;
yterm := yterm * xj;
END;
END;
{ determine the regression coefficients }
FOR j := uno TO nterms DO
FOR k := uno TO nterms DO
BEGIN
n := j + k - uno;
arry[j,k] := sumx[n];
END;
delta := determ(arry,nterms);
IF delta = zero THEN
bad_data(1)
ELSE
BEGIN
FOR l := uno TO nterms DO
BEGIN
FOR j := uno TO nterms DO
BEGIN
FOR k := uno TO nterms DO
BEGIN
n := j + k - uno;
arry[j,k] := sumx[n];
END;
arry[j,l] := sumy[j];
END;
a^[l] := determ(arry,nterms) / delta
END;
{ calculate r and standard error of estimate }
srs := zero;
sum_y := zero;
sum_y2 := zero;
FOR j := uno TO npts DO
BEGIN
yj := y^[j];
xj := x^[j];
ycalc := zero;
x_power := one;
FOR k := uno TO nterms DO
BEGIN
ycalc := ycalc + a^[k] * x_power;
x_power := x_power * xj;
END;
residual := ycalc - yj;
srs := srs + SQR(residual);
sum_y := sum_y + yj;
sum_y2 := sum_y2 + (yj * yj);
END;
sigmay2 := (sum_y2 - SQR(sum_y)/npts)/(npts-1);
IF sigmay2 <= zero THEN
bad_data(2)
ELSE
BEGIN
se2 := srs/(npts - 2);
dummy := se2/sigmay2;
IF dummy <= one THEN
BEGIN
r := SQRT(one - dummy);
IF (nterms = 2) AND (a^[nterms] < zero) THEN
r := -r;
se := SQRT(se2);
END
ELSE
bad_data(3);
END;
END;
END;
PROCEDURE smooth121(n:WORD; VAR y: single_array_pointer);
{ Author : Norton Associates
Purpose: Smooth data via binomial 121
Version: 1.0
Date : 5 May 1990 }
VAR
nmax,j : WORD;
y_old,
new_y : SINGLE;
BEGIN
{ check to see if there is enough data }
IF n < 3 THEN
BEGIN
WRITELN('smooth121> no sample data');
EXIT;
END;
nmax := n - 1;
y_old := y^[1];
{ smooth data }
FOR j := 1 TO nmax DO
BEGIN
new_y := (y_old + c2 * y^[j] + y^[j+1]) * i_4;
y_old := y^[j];
y^[j] := new_y;
END;
{ edge effect }
y^[n] := (y_old + 3.0 * y^[n])* i_4;
END;
PROCEDURE smooth14641(n:WORD; VAR y: single_array_pointer);
{ Author : Norton Associates
Purpose: Smooth data via binomial 14641
Version: 1.0
Date : 5 May 1990 }
VAR
nmax,j : WORD;
y_2,y_1,
new_y : SINGLE;
BEGIN
{ see if there is enough data }
IF n < 5 THEN
BEGIN
WRITELN('smooth14641> no sample data');
EXIT;
END;
nmax := n - 2;
y_2 := y^[1];
y_1 := y^[1];
FOR j := 1 TO nmax DO
BEGIN
new_y := (y_2 + c4*y_1 + c6*y^[j] + c4*y^[j+1] + y^[j+2])* i_16;
y_2 := y_1;
y_1 := y^[j];
y^[j] := new_y;
END;
new_y := y^[n-1];
{ edge effects }
y^[n-1] := (y_2 + c4*y_1* c6*y^[n-1] + c5*y^[n]) * i_16;
y^[n] := (y_1 + c4*new_y + c11*y^[n]) * i_16;
END;
PROCEDURE smoothcurve(n:WORD; VAR y: single_array_pointer);
{ Author : Norton Associates
Purpose: Smooth curvelinear data
Version: 1.0
Date : 5 May 1990 }
VAR
nmax,j : WORD;
y_2,y_1,
new_y : SINGLE;
BEGIN
{ check to see if there is enough data }
IF n < 6 THEN
BEGIN
WRITELN('smoothcurve> no sample data');
EXIT;
END;
nmax := n - 2;
y_2 := y^[1];
y_1 := y^[1];
FOR j := 1 TO nmax DO
BEGIN
new_y := ( c3*(y_2 + y^[j+2]) + c12*(y_1 + y^[j+1]) + c17*y^[j])* i_35;
y_2 := y_1;
y_1 := y^[j];
y^[j] := new_y;
END;
new_y := y^[n-1];
{ edge effects }
y^[n-1] := (y_2 + c4*y_1 * c6*y^[n-1] + c5*y^[n]) * i_16;
y^[n] := (y_1 + c4*new_y + c11*y^[n]) * i_16;
END;
PROCEDURE corcoef(n: WORD; x,y:single_array_pointer; VAR r:SINGLE);
{ Author : Norton Associates
Purpose: Calculate correlation coefficient based upon
product moment algrothim
Version: 1.0
Date : 5 May 1990 }
VAR
sumx,sumy,sumxy,sumy2,sumx2 : DOUBLE;
j : WORD;
BEGIN
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('corcoef> no sample data');
EXIT;
END;
{ zero out some data }
sumx := 0.0;
sumy := 0.0;
sumxy := 0.0;
sumy2 := 0.0;
sumx2 := 0.0;
{ calculate sum of squares }
FOR j := 1 TO n DO
BEGIN
sumx := sumx + x^[j];
sumy := sumy + y^[j];
sumxy := sumxy + x^[j]*y^[j];
sumx2 := sumx2 + SQR(x^[j]);
sumy2 := sumy2 + SQR(y^[j]);
END;
{calculate correlation coefficient }
r := ( (n*sumxy) - (sumx*sumy) )/
( ( SQRT((n*sumx2 - SQR(sumx))) * (SQRT(n*sumy2 - SQR(sumy))) ));
END;
PROCEDURE mulreg(n : WORD;
y,x,z : single_array_pointer;
VAR a : single_array_pointer;
VAR r : SINGLE;
VAR se: SINGLE);
{ Author : Norton Associates
Purpose: Determine via least squares the regresion
and correlation coefficients for multiple regression
Version: 1.0
Date : 5 May 1990 }
VAR
smx,smy,smz,smz2,smx2,smxz,smxy,smyz,
delta,srs,sum_y,sum_y2,
ycalc,x_power,residual : EXTENDED;
sumy : ARRAY[1..maxorder] OF EXTENDED;
arry : arry_type;
j,nmax,k,l : WORD;
xj,yj,zj : SINGLE;
sigmay2,se2,dummy : EXTENDED;
PROCEDURE bad_mul_data(nn : WORD);
VAR
j : WORD;
BEGIN
r := zero;
se := zero;
FOR j := 1 TO 3 DO
a^[j] := zero;
WRITELN('mulreg> Determinate or correlation = 0.0 ',nn);
END;
PROCEDURE coeff(VAR coef : SINGLE; kk : WORD ; tarry : arry_type);
VAR
i : WORD;
BEGIN
FOR i := 1 TO 3 DO
BEGIN
tarry[i,kk] := sumy[i];
IF kk > 1 THEN
tarry[i,kk - 1] := arry[i,kk - 1];
END;
coef := determ(tarry,3) / delta;
END;
BEGIN
{ check to see if there is enough data }
IF n <= 3 THEN
BEGIN
WRITELN('mulreg> sample data less than 4');
EXIT;
END;
{ zero data out }
r := zero;
FOR j := 1 TO 3 DO
a^[j] := zero;
smy := zero;
smz := zero;
smx := zero;
smx2 := zero;
smz2 := zero;
smxz := zero;
smyz := zero;
smxy := zero;
{ calculate sum of squares for data }
FOR j := 1 TO n DO
BEGIN
yj := y^[j];
xj := x^[j];
zj := z^[j];
smy := smy + yj;
smx := smx + xj;
smz := smz + zj;
smx2 := smx2 + SQR(xj);
smz2 := smz2 + SQR(zj);
smxz := smxz + (xj*zj);
smxy := smxy + (xj*yj);
smyz := smyz + (yj*zj);
END;
{ set up linear simultaneous equations}
sumy[1] := smy;
sumy[2] := smxy;
sumy[3] := smyz;
arry[1,1] := n;
arry[2,1] := smx;
arry[3,1] := smz;
arry[1,2] := smx;
arry[2,2] := smx2;
arry[3,2] := smxz;
arry[1,3] := smz;
arry[2,3] := smxz;
arry[3,3] := smz2;
delta := determ(arry,3);
IF delta = zero THEN
bad_mul_data(1)
ELSE
BEGIN
coeff(a^[1],1,arry);
coeff(a^[2],2,arry);
coeff(a^[3],3,arry);
END;
{ calculate r and standard error of estimate }
srs := zero;
sum_y := zero;
sum_y2 := zero;
FOR j := 1 TO n DO
BEGIN
yj := y^[j];
xj := x^[j];
zj := z^[j];
ycalc := a^[1] + a^[2]*xj + a^[3]*zj;
residual := ycalc - yj;
srs := srs + SQR(residual);
sum_y := sum_y + yj;
sum_y2 := sum_y2 + SQR(yj);
END;
sigmay2 := (sum_y2 - SQR(sum_y)/n)/(n - 1);
IF sigmay2 <= zero THEN
bad_mul_data(2)
ELSE
BEGIN
se2 := srs/(n - 3);
dummy := se2/sigmay2;
IF dummy <= one THEN
BEGIN
r := SQRT(one - dummy);
se := SQRT(se2);
END
ELSE
bad_mul_data(3);
END;
END;
PROCEDURE autocor(n:WORD; x:single_array_pointer; lag :WORD; VAR auto : single_array_pointer);
{ Author : Norton Associates
Purpose: Determine correlation coefficients for auto lagging of the data
Version: 1.0
Date : 5 May 1990 }
VAR
k,j,
num : WORD;
y : single_array_pointer;
r : SINGLE;
BEGIN
{ check to see if there is enough data }
IF n < dos THEN
BEGIN
WRITELN('autocor> no sample data');
EXIT;
END;
{ lag value can not be less than the number of data points }
IF lag > n THEN
BEGIN
WRITELN('autocor> lag can not exceed number in sample data');
EXIT;
END;
create_single_array(n,y);
FOR j := 1 TO lag DO
BEGIN
num := n - j + 1;
MOVE(x^[j],y^[1],(num*SIZEOF(x^[1])));
corcoef(num,x,y,r);
auto^[j] := r;
END;
delete_single_array(n,y);
END;
FUNCTION movavg(n: WORD; a : single_array_pointer; ma:WORD; k : WORD) : SINGLE;
{ Author : Norton Associates
Purpose: Determine the moving average of data for a selected data set
Version: 1.0
Date : 5 May 1990 }
VAR
dum : SINGLE;
j : WORD;
BEGIN
{ check to see if index is out of range }
IF (k > n) OR ( k < 1) THEN
BEGIN
WRITELN('movavg> sorry index is out of range');
EXIT;
END;
{ check to see if the index is less than the lag number }
IF k < ma THEN
BEGIN
WRITELN('movavg> number of points less than moving average: movavg = 0');
movavg := 0.0;
END
ELSE
BEGIN
dum := 0.0;
FOR j := (k - ma + uno) TO k DO
dum := dum + a^[j];
movavg := dum / ma ;
END;
END;
BEGIN
END.